## [1] <NA> 1 2 3 4
## Levels: 1 2 3 4
On récupère les données suivantes :
obsdata, segdateet predatadistdata, une jointure entre predata et obsdata sur le segmentOn centre et réduit les covariables (présentes dans distdata, segdata et predata).
## CHL_4w_mea CHL_4w_sd SST_4w_mea SST_4w_sd POC_4w_mea depth distCoast
## mean 1.3016453 0.4070661 15.865554 0.8412628 211.17688 51.14547 38510.47
## sd 0.8536337 0.3629038 3.251786 0.4386712 94.86176 36.32341 30059.10
## dist200 slopeP
## mean 88112.13 28913.18
## sd 32477.64 27947.08
On recalcule X et Y en lambert 93 à cause d’un souci dans les données.
predata
distdata
gridata
## # A tibble: 4 x 2
## # Groups: session [4]
## session dates_session
## <fct> <chr>
## 1 1 2019-02-12 + 2019-02-13 + 2019-02-14 + 2019-02-15 + 2019-02-26 + 2019~
## 2 2 2019-05-30 + 2019-05-31 + 2019-06-01 + 2019-06-02 + 2019-06-03
## 3 3 2019-07-31 + 2019-08-01 + 2019-08-02 + 2019-08-05 + 2019-08-07 + 2019~
## 4 4 2019-10-25 + 2019-10-26 + 2019-10-27 + 2019-10-28 + 2019-11-17 + 2019~
On cherche à savoir quelle est la meilleure fonction de détection.
## [1] "Itération 1 sur 8 terminée."
## [1] "Itération 2 sur 8 terminée."
## [1] "Itération 3 sur 8 terminée."
## [1] "Itération 4 sur 8 terminée."
## Error :
## gosolnp-->Could not find a feasible starting point...exiting
##
## [1] "Itération 5 sur 8 terminée."
## [1] "Itération 6 sur 8 terminée."
## [1] "Itération 7 sur 8 terminée."
## [1] "Itération 8 sur 8 terminée."
## observerId seaState formula key AIC
## 1 0 0 ~1 hn -293.4346
## 2 0 1 ~seaState hn -290.1428
## 3 1 0 ~observerId hn -288.9999
## 4 1 1 ~observerId+seaState hn -288.8973
## 5 0 0 ~1 hr -297.3398
## 6 0 1 ~seaState hr -297.3362
## 7 1 0 ~observerId hr -295.5575
## 8 1 1 ~observerId+seaState hr -295.8082
On choisit la fonction de détection avec l’AIC le plus faible. Ici, il s’agit de la fonction avec la formule ~1 et la key hr. A cause d’une erreur pour ajuster cette fonction seulement, on préferera prendre la fonction avec la formule ~seaState et la key hr, qui a un AIC très très proche. On sélectionne donc la meilleure fonction de détection : detfc.sea.hr
##
## Summary for distance analysis
## Number of observations : 96
## Distance range : 0 - 0.282
##
## Model : Hazard-rate key function
## AIC : -297.3362
##
## Detection function parameters
## Scale coefficient(s):
## estimate se
## (Intercept) -2.2314812 0.1847034
## seaState 0.1731484 0.1239956
##
## Shape coefficient(s):
## estimate se
## (Intercept) 1.473509 0.242613
##
## Estimate SE CV
## Average p 0.5468238 0.04173954 0.07633087
## N in covered region 175.5593056 18.07215909 0.10294048
On calcule manuellement les courbes de détection en fonction de l’état de la mer.
# On récupère l'ensemble des paramètres du modèle d'abord
detfc_par <- detfc.sea.hr$ddf$par
# scale.value : fonction de mise à l'échelle des paramètres du modèle (matrice et exponentielle)
scale.value <- function (param, z){
exp(as.matrix(z) %*% param) }
# key.fct.hz : On entre une distance, un sigma (scale), un b (shape) et on obtient une probabilité de détection pour une hazard-rate
key.fct.hr <- function (distance, sigma, b){
return(1 - exp(-(distance / sigma) ^ (-b))) }
# On créé un vecteur des distances de 0 à 0.3 km
dist <- seq(from = 0, to = 0.3, length.out = 96)
# On récupère le paramètre shape (b) et on le met à l'échelle : 100 fois exp(b)
b <- scale.value(detfc_par["V1"], matrix(1, nrow = 96, 1))
# On va ensuite récupérer le paramètre scale (sigma)
# C'est plus compliqué ici, car sigma varie selon les covariables de détection, donc seaState ici.
# sigma.0 : intercept
sigma0 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState = 0
z = matrix(c(rep(1, 96), rep(0, 96)), ncol = 2))
# sigma.1 : intercept + seaState toujours à 1
sigma1 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState = 1
z = matrix(rep(1, 96*2), ncol = 2))
# sigma.2 : intercept + seaState toujours à 2
sigma2 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState constant à 2
z = matrix(c(rep(1, 96), rep(2, 96)), ncol = 2))
# sigma.3 : intercept + seaState toujours à 3
sigma3 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState constant à 2
z = matrix(c(rep(1, 96), rep(3, 96)), ncol = 2))
# A partir de tous ces paramètres et du vecteur des distances, on calcule la probabilité de détection
proba_detec0 <- key.fct.hr(distance = dist, sigma = sigma0, b = b)
proba_detec1 <- key.fct.hr(distance = dist, sigma = sigma1, b = b)
proba_detec2 <- key.fct.hr(distance = dist, sigma = sigma2, b = b)
proba_detec3 <- key.fct.hr(distance = dist, sigma = sigma3, b = b)On va ensuite générer le graphique correspondant.
On enregistre ce graphique dans le dossier img en 2 versions, en anglais et en français, pour le poster.## png
## 2
## png
## 2
Nous avons aussi tenté en mettant les modalités beaufort 0-1 ensemble et 2-3 ensemble, pour voir si nous avions toujours une meilleure détection avec une mer plus agitée.
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 25992 individuals, described by 9 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"
On utilise la fonction selec_dsm_aic_fwd qui va permettre de sélectionner de manière forward les covariables à inclure dans le modèle dsm.
Cette fonction prend en arguments segdata, obsdata, la fonction de détection et un vecteur de toutes les covariables à tester. C’est une fonction récursive, elle peut donc aussi prendre en argument un vecteur des covariables déjà sélectionnées.
On sélectionne les covariables sur l’ensemble des données, sans distinction de session. Ici, les AIC sont données pour
availability = 1. L’algorithme pour les autres valeurs de disponibilité sélectionne les même covariables, seul l’AIC change.
## [1] "----- Sélection de la covariable numéro 1 -----"
## [1] "----- Sélection de la covariable numéro 2 -----"
## [1] "----- Sélection de la covariable numéro 3 -----"
## [1] "----- Sélection de la covariable numéro 4 -----"
## [1] "L'AIC ne diminue plus en ajoutant une 4ème covariable"
## $dsm.selec
##
## Density surface model
## Response : count
##
## Detection function : Hazard-rate key function
##
## Formula: count ~ s(SST_4w_mea) + s(X, Y) + s(CHL_4w_mea) + offset(off.set)
##
## Fitting engine: gam
##
##
## $aic
## [1] 1757.999
##
## $formule.dsm
## [1] "count ~ s(SST_4w_mea) + s(X, Y) + s(POC_4w_mea)"
##
## $suivi
## iter formule aic
## 1 1 count ~ s(SST_4w_mea) 1842.067
## 2 2 count ~ s(SST_4w_mea) + s(X, Y) 1763.128
## 3 3 count ~ s(SST_4w_mea) + s(X, Y) + s(CHL_4w_mea) 1757.999
dsm pour \(availability\) dépendante de on-shelf et off-shelf : On note “on-shelf” quand la profondeur est inférieure à 150m, et “off-shelf” si la profondeur est supérieure à 150m.
\[availability_{off-shelf}=0,1357617\] \[availability_{on-shelf}=0,6332016\]
On créé deux tableaux de données predata_tmp dont les valeurs correspondent aux valeurs de predata pour chaque session.
## Session Disponibilité Min Estimation Max
## 1: 2 1 2948.3582 4441.895 5935.433
## 2: 3 1 348.7289 1536.801 2724.874
## 3: 2 0.41 9170.5957 10664.133 12157.670
## 4: 3 0.41 2515.4844 3703.557 4891.629
## 5: 2 on-shelf/off-shelf 5452.3439 6945.881 8439.418
## 6: 3 on-shelf/off-shelf 1194.0925 2382.165 3570.237
## Session Disponibilité Estimation Plus.ou.moins
## 1: 2 1 4441.895 1493.537
## 2: 3 1 1536.801 1188.072
## 3: 2 0.41 10664.133 4266.205
## 4: 3 0.41 3703.557 3566.324
## 5: 2 on-shelf/off-shelf 6945.881 2531.880
## 6: 3 on-shelf/off-shelf 2382.165 2074.952
On créé maintenant toutes les cartes statiques. Ce sont celles que l’on retrouve sur le poster. Elles vont également être enregistrées dans le dossier img.
## png
## 2
## png
## 2
## png
## 2
## png
## 2
## png
## 2
## png
## 2
On obtient aussi une carte interactive par session, qui sont celles que l’on retrouve dans l’application.
On enregistre plusieurs images dans le dossier img : par session et biais de disponibilité, on a une image avec toutes les splines et une image avec X et Y en plus grand.
Voici par exemple l’image pour la sessio 2 (printemps), un biais de 1 et toutes les covariables.
On récupère les predata avec les covariables non centrées-réduites.
On retrouve toutes les cartes ci-dessous :
On exporte certains résultats en Rdata qui seront utilisés pour l’application shiny.